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Abstract 

We present an improved method for topology optimization with both adaptive mesh refinement and 
derefinement. Since the total volume fraction in topology optimization is usually modest, after a few 
initial iterations the domain of computation is largely void. Hence, it is inefficient to have many small 
elements, in such regions, that contribute significantly to the overall computational cost but contribute 
little to the accuracy of computation and design. At the same time, we want high spatial resolution 
for accurate three-dimensional designs to avoid postprocessing or interpretation as much as possible. 
Dynamic adaptive mesh refinement (AMR) offers the possibility to balance these two requirements. We 
discuss requirements on AMR for topology optimization and the algorithmic features to implement them. 
The numerical design problems demonstrate (1) that our AMR strategy for topology optimization leads 
to designs that are equivalent to optimal designs on uniform meshes, (2) how AMR strategies that do 
not satisfy the postulated requirements may lead to suboptimal designs, and (3) that our AMR strategy 
significantly reduces the time to compute optimal designs. 

keywords: adaptive mesh refinement, topology optimization, iterative solvers. 

1 Introduction 

Topology optimization is a powerful structural optimization method that combines a numerical solution 
method, usually the finite element method, with an optimization algorithm to find the optimal material 
distribution inside a given domain (THl [23 HSl IH] • In designing the topology of a structure we determine 
which points in the domain should be material and which points should be void. However, it is known that, 
in the continuum setting, topology optimization leads to designs with intermediate densities. So continuous 
values between and 1 replace discrete values (0 or 1) to represent the relative densities, and some form of 
penalization is used to obtain designs with almost discrete 0/1 material density distribution [3]. 

In topology optimization, problems are solved most commonly on fixed uniform meshes with a relatively 
large number of elements in order to achieve accurate designs [T71 [T^]. However, as void and solid regions 
appear in the design, it is more efficient to represent the holes with fewer large elements and the solid regions, 
especially the material surface, with more fine elements. Since the shape and position of holes and solid 
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regions are initially unknown, the most economical mesh representation for the design is unknown a priori. 
Therefore, adaptive mesh refinement (AMR) is very suitable for topology optimization. The purpose of AMR 
for topology optimization is to get the design that would be obtained on a uniformly fine mesh, but at a much 
lower computational cost by reducing the total number of elements and having fine elements only where and 
when necessary. 

Highly accurate designs on uniform meshes may require so many elements that the solution of the 
optimization problem becomes intractable. However, AMR leads to high resolution in the mesh only when 
and where necessary. This makes it possible to obtain accurate designs with a modest number of elements 
and hence with a reasonable cost. Even when a design on a uniform mesh is computationally feasible, AMR 
tends to reduce the computational cost by reducing the required number of elements and by improving the 
conditioning of linear systems arising from the finite element discretization. Obviously, we do not want the 
use of AMR or the AMR procedure to alter the computed designs. However, there is a risk of this, since 
the mesh influences the computed deformations and sensitivities. It is imperative then that the solutions 
from the finite element analysis using AMR must be as accurate as those obtained on a uniform fine mesh. 
Moreover, the final design must be governed by accurate sensitivities corresponding to those obtained 
on the finest mesh. If coarse mesh solutions drive or limit the design, suboptimal designs may result when 
designs optimal on a coarser mesh differ substantially from the optimal design on a (much) finer mesh. We 
will demonstrate that this occurs in quite simple cases. The early work in this area, though leading to 
acceptable designs in specific instances, does not satisfy these properties. We will propose relatively simple 
but essential changes to these methodologies that lead to AMR-based designs that are equivalent (up to 
some small tolerance) to designs on uniform fine meshes. In addition, our approach leads (in principle) to an 
overall more efficient method as we reduce the total number of elements further. The topology optimization 
may lead to a sequence of (intermediate) structures requiring high mesh resolution in different parts of the 
computational domain. Therefore, it is important to (1) allow the meshes at all levels to change continually 
(dynamic) and (2) to allow both mesh refinement and derefinement [23j . Derefinement is important for 
efficiency when the initial discretization needs to include relatively small elements in certain regions. This 
is important in a number of cases, which are elaborated upon below. 

In the next section, we provide an assessment of previous AMR strategies, namely the implementations by 
Costa and Alves [8] and Stainko [22]. In Section |3l we provide a brief introduction to topology optimization. 
Next, in Sectional we state the purpose of our AMR strategy for topology optimization and explain the 
requirements it poses. Based on these requirements, we propose a more robust and dynamic AMR strategy. 
We describe a number of implementation issues of our AMR strategy in Section [Bj We briefiy discuss 
the iterative solution of the large sparse linear systems arising in the finite element analysis in Section [5] 
In Section [71 we show numerical experiments that demonstrate the robustness and efficiency of our AMR 
strategy. The first experiment also explains why the refinement strategies by Costa and Alves [8 and Stainko 
[22] may lead to suboptimal designs. Finally, in Section [8] we present conclusions about our AMR strategy 
for topology optimization algorithms, and we mention some directions for future work. 

2 Assessment of Previous AMR Strategies 

Little research has been done in applying AMR to topology optimization. So, we start by briefly discussing 
two recent, important, papers in this area. The AMR method by Costa and Alves [8] goes through a prede- 
termined, fixed sequence of optimizations and subsequent mesh refinements (they do not use derefinements), 

^For comparison everywhere in this paper the element size for the uniform fine mesh is the same as the element size at the 
highest level of refinement in the AMR, mesh. 
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always using (or assuming) a converged solution on a 'coarse mesh' to guide the refinement of that mesh 
and start the optimization on the next 'fine mesh'. Coarse meshes and the solutions on these coarse meshes 
are never revisited or updated after generating the next finer mesh. The method aims at refining the coarse 
mesh design. Hence the region with the fine(r) mesh that contains the material boundary is always confined 
to the coarse mesh region that has been determined before using only coarse mesh calculations. After a fixed 
number of optimization steps on a given mesh, they refine all material elements (density 0.5 or larger) and 
elements on the boundary between material elements and void elements (density less than 0.5). Furthermore, 
they refine elements that do not satisfy certain quality or solution error criteria. In addition, there are a 
few special cases that also lead to refinement. These refinements lead to accurate finite element solutions in 
material regions, a high mesh resolution on the material boundary and, therefore, accurate representation 
of this boundary, and larger elements in void regions reducing the amount of work. However, as reported 
by the authors, the 'optimal design' found by the method depends on the initial mesh and is not the same 
as the optimal design found using a uniform fine mesh [F. Although the authors do not report this, we 
conjecture that the design found using the adaptively refined mesh is not an 'optimal design' on the uniform 
mesh, that is, it has higher compliance than the solution obtained on the uniform fine mesh. See also our 
numerical experiments below. Finally, only two-dimensional designs are treated, but conceptually we expect 
their algorithm to work similarly in three-dimensional designs. 

Stainko follows a slightly different approach with respect to the refinements [22] . Mesh refinement is done 
only along the material boundary as indicated by the (regularization) filter. So, elements completely inside 
a material region or a void region are not refined. In principle this leads to a smaller number of elements 
and hence a reduced computation time. However, Stainko's procedure also progresses strictly from coarser 
meshes to finer meshes, and a coarse mesh is never updated after mesh refinement. So, just as in the finest 
mesh, which contains the material boundary, is always confined to regions defined by (all) earlier refinements 
(all refinements are nested), each of which is based only on the corresponding coarse mesh computations. 
Stainko does not test whether the designs obtained are the same as those obtained on uniformly fine meshes; 
however, our experiments below show that, again, the designs will be depend on the initial mesh (resolution) 
and are not the same as optimal designs on the uniformly fine mesh at the maximum refinement level. 

These approaches share two important choices that may lead to problems. First, both approaches solve 
the design problem on a fixed mesh until convergence before carrying out mesh refinement. After refinement 
on a given level, the mesh on that level remains fixed for the remainder of the optimization, and all further 
refinements are therefore constrained by the converged coarser level solutions. This works well in terms of 
refining the design, but for many design problems the optimal solution on a uniform fine(st) mesh is quite 
different from the converged solution on a coarser mesh. In that case, mesh refinement based only on the 
coarser level solution will erroneously confine the solution on the finer mesh to a smooth version of the 
coarser level solution. Therefore, the approaches proposed in [22l |8] may lead to suboptimal designs, as we 
will show in our numerical experiments. 

Second, both approaches use only refinement but no derefinement, which may lead to inefficiencies. First, 
for designs with thin structures, the initial, coarsest, mesh must be fine enough to give a reasonable result. If 
fine elements that are no longer required cannot be removed as the design evolves, then more computational 
work than necessary will be performed. Second, in topology optimization approaches that use filtering for 
regularization, for an accurate design requiring a high resolution mesh, the (appropriate) filter will not work 
on the coarser meshes, because the filter radius, which should be a physical feature size independent of the 
mesh, will typically be too small. Hence, we must start with a relatively fine mesh. However, after a modest 
number of optimization steps, large regions will likely have become void and fine elements could be removed 
without problems. Again substantial computational overhead results from having to work with too many 
fine elements. Third, any AMR strategy that allows changes in the design beyond previously computed 



3 




(Final iteration) 

Figure 1 : Overview of the topology optimization algorithm with dynamic AMR. 

coarse level designs and refines the mesh to accommodate such changes will be inefficient if fine elements in 
void regions cannot be removed. 

Therefore, a more robust and efficient refinement strategy is needed. Hence, we propose a dynamic 
meshing strategy that includes both mesh refinement and derefinement everywhere in the computational 
domain. Our improved AMR strategy has two main components. First, we extend the refinement criteria 
from [8 , refining all material elements and elements on the boundary, but with an additional layer of 
refinements around the material boundary (in the void- region) . The thickness of the layer is a parameter. 
This way the fine level design can change shape arbitrarily in optimization steps between mesh refinements. 
Second, our AMR method updates coarse and fine meshes continually, so that small changes arising in the 
more accurate computations on finer meshes can change the shape of the design arbitrarily in the course 
of the optimization; the fine(r) meshes move with the material boundary. This means that our designs are 
really based on the accurate fine mesh computations and are not confined to regions fixed by earlier coarse 
mesh computations. Since we do continual mesh adaptation, we may have fine elements in regions that 
have become void at some point. Derefinement will remove those fine elements. Further details are given in 
the next subsection. This approach leads to designs on AMR meshes at greatly reduced cost that are the 
same as the designs that would have been obtained on a uniform fine mesh of the highest resolution (within 
a small tolerance) . We will demonstrate this experimentally in section [T] Our approach also allows us to 
start with coarser meshes, since the coarse mesh solution does not need to be a good approximation to the 
final solution. Even when we start with a finer mesh for faster convergence or proper functioning of the 
regularization filter, derefinement allows us to remove fine elements that have become void (which tends to 
happen quite quickly). 
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3 A Brief Topology Optimization Review 



In topology optimization we solve for the material distribution in a given computational design domain f2. 
The topology optimization problem we consider here is to minimize the compliance of a structure under 
given loads as a function of the material distribution. To solve this problem numerically, we discretize the 
computational domain using finite elements, where we usually use a lower order interpolation for the density 
field (material distribution) than for the displacement field. The most common approach (also employed 
for this paper) is to use (bi-,tri-)linear interpolation for the displacement field and constant density in each 
element. The compliance minimization problem after finite element discretization is defined as 

min f^u (1) 

PeG[ft,ll,Ve 

{K{p)u = f for a; £ O \ Oq, 
u = Uq for X G r^Oj 

EePeVe<Vo, 

where pe is the density in element e, p is the vector of element densities, K is the stiffness matrix, a function 
of the discretized density field (p), is the volume of element e, Vq is a maximum volume (fraction) allowed 
for the design, and f2o is the part of the domain where the displacement is prescribed. To avoid singularity 
of the stiffness matrix, we enforce a small positive lower bound on the element density, typically 10"'^. 

As mentioned in the introduction, our discrete model must drive the continuous material distribution as 
much as possible to a 0/1 solution. We use the Solid Isotropic Material with Penalization (SIMP) method 
to make the undesirable intermediate densities between /9 (replacing 0) and 1 unfavorable [6 . In this case, 
the elasticity tensor is defined as a function of the element density, 

Ee - P^eEo, (2) 

where p is the penalization parameter. With p > 1, intermediate densities are unfavorable as they provide 
relatively little stiffness compared with their material cost. A common choice is p = 3, which results in 
intermediate material properties that satisfy the Hashin-Shtrikman bound for composite materials [lOj . To 
avoid problems with local minima, we usually apply continuation on the parameter p, that is, we start with 
p — 1 and slowly increase p as the design converges. 

The general scheme for topology optimization using AMR is illustrated in Figure [TJ First, we set up the 
geometry, the finite element (FE) mesh, the loading and boundary conditions, and we initialize the density 
distribution p. Then, we start the optimization loop. In this loop, we assemble and solve the equilibrium 
equations K{p)u = / in ([T]) using the FE discretization and a linear solver. Next, in the sensitivity analysis, 
we compute the derivatives of the objective function with respective to the design variables, dc/dpe- After 
this, we can apply an optional low-pass filter to remedy the checkerboard problem |18 [ I19 | [2T|. which can be 
also addressed by an alternative minimum length scale approach [5] . In the next step, we compute an update 
of the design variables. There are various optimization algorithms applicable to topology optimization. For 
this paper, we use Optimality Criteria (OC), a simple approach based on a set of intuitive criteria [6l |4]. 
After updating the design variables using a chosen optimization algorithm, we check the convergence of the 
design. Under certain conditions, to be discussed next, dynamic mesh adaptation is carried out before the 
(next) finite element analysis. 
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4 A Dynamic AMR Strategy 



We base our algorithmic choices on a set of requirements on AMR codes for topology optimization. As 
stated above, the purpose of AMR for topology optimization is to get the design that would be obtained on 
a uniform fine mesh, but at a much lower computational cost by reducing the total number of elements and 
having fine elements only where (and when) necessary. 

First, since the finite element analysis and the computation of sensitivities drive the changes in material 
distribution, they should be as accurate as on the uniform fine mesh. Therefore, we need a fine mesh that 
covers at least the material region and the boundary. Since the void regions have negligible stiffness they do 
not infiuence the (intermediate) linear finite element solutions and resulting sensitivity computations. Thus 
we do not need a fine mesh inside the void region, and we can use a refinement criterion similar to that of 
Costa and Alves [8 . At this point we focus on refinement and dcrcfinemcnt for shape only. Therefore, we are 
conservative with respect to accuracy, and we expect that, in future implementations, good error indicators 
will lead to further efficiency gains, in particular because of derefinement in solid material regions. 

Second, the accurate computations on the finest level should drive the changes in the material distribution. 
This requires continual mesh adaptation so that computational results after refinements can drive updates 
to the material distribution, and designs are not confined by earlier coarse grid results. This also means that 
as the material region moves close to the boundary between fine and coarse(r) mesh, additional refinements 
allow for further evolution. 

Third, we need to ensure that the design can change sufficiently in between mesh updates. Therefore, 
we maintain a layer of additional refinements around the material region (in the void region) and carry out 
continual mesh adaptation. Due to the additional layer of refinements and continual mesh updates, the design 
can change arbitrarily following the fine grid computations and resulting sensitivities, and it is not confined 
by earlier coarse grid results. To ensure that the design accurately reflects the fine mesh computations, we 
allow rapid refinements of the mesh early on when voids and material regions (and hence the boundary) 
develop, rather than delay refinements until later stages, when a suboptimal design might have developed. 

Fourth, since the design can change substantially from its estimate on a coarse mesh, we may have fine 
elements in void regions. Those elements must be removed for efficiency, requiring derefinements of the mesh. 
To facilitate our strategy of continual mesh refinement and derefinement, we use a hierarchical representation 
of adaptive meshes. 

We will now state our refinement and derefinement strategy in more detail. We adapt the mesh when 
CASE (i): 

1. the relative change in the compliance is smaller than a given threshold, and 

2. a given minimum number of optimization steps have been carried out since the last mesh update, 

or when 
CASE (ii): 

3. a given maximum number of optimization steps have been carried out without meeting conditions 1 
and 2. 

Condition 1 corresponds to a common convergence criterion for topology optimization that the maximum 
change in the design variables is smaller than a certain tolerance (which we usually set to 0.01). This 
condition is satisfied when the solution is near a local minimum, which might be caused by a no longer 
appropriate mesh. In that case, we must adapt the mesh to allow the design to change further. If the local 
minimum is not an artifact of the mesh, the design will remain the same after mesh adaptation. Condition 2 
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Figure 2: Refinement criteria for void element. Element a is marked for refinement, because it has solid 
elements within distance Tamr! element h is marked for derefinement. 

prevents refinement and derefinement fi'om happening too fi'equently. This is important as the solution needs 
to adapt to the changed mesh, so that the computed sensitivities refiect the design and not mesh artifacts. 
This also limits the cost of mesh adaptation. In our experiments this minimum number of optimization steps 
is set to five based on experience, and other values are possible. Regarding condition 3, in our experiments 
we adapt the mesh at least every ten optimization steps. This condition leads to faster convergence because 
it ensures that the mesh is appropriate for the material distribution. Using these conditions, we can start 
with a fairly coarse mesh, and we may carry out mesh (de)refinement before the design converges on any 
mesh if necessary. 

We adapt our mesh according to the following procedure. 

1. Mark all the elements for refinement or derefinement based on the following criteria: 

• If element e is solid, i.e., G [ps, 1], where ps is a chosen density threshold, or element e is within 
a given distance Tamr from a solid element we mark it for refinement. 

• If element e is void, i.e., pe £ [p^^ps], and there are no solid elements within distance Tamr, we 
mark element e for derefinement. See Figure [5J 

2. Check compatibility for the mesh that will be generated and make the following adjustments in two 
sweeps over all elements: 

• In the first sweep, we unmark elements marked for derefinement, if they have a sibling (an element 
generated by the same refinement) that is not marked for derefinement. 

• In the second sweep, we unmark elements marked for derefinement, if derefinement would lead 
to level two or higher edge incompatibility. We allow only level one incompatibility; see Figure |3] 
and the discussion in Section [S) 

The above refinement criteria result in a layer of fine elements on the void side of the solid/void interface 
that allows the material to be redistributed locally. If a material boundary moves near the fine/coarse 
element interface, mesh refinement creates a new layer of fine elements around the current material surface 
to allow further local redistribution of the material. On the other hand, if some fine elements become void, 
these fine elements are removed by derefinement to keep the optimization efiicient. 
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5 Iterative Solution Scheme 



Although AMR significantly reduces the number of DOFs in the finite element simulation, we still have to 
solve a sequence of large linear systems, especially for three-dimensional designs. Moreover, because of the 
large difference in density between distinct parts of the computational domain and the elasticity tensor given 
by (21), with p ^ 3 toward the end of the optimization, the linear systems are very ill-conditioned. Hence, 
proper preconditioning is essential. In |24j . we showed how to precondition the linear systems arising in 
topology optimization, and we also used Krylov subspace recycling il3j to reduce the number of iterations 
over multiple linear systems. We briefly mention the main ideas here. 

To remedy the serious ill-conditioning in topology optimization problems, we explicitly rescale each 
stiffness matrix such that the diagonal coefficients are all equal, as is the case for a problem with homogeneous 
density. We rescale the stiffness matrix K by multiplying with a diagonal matrix on both sides, 

where D is the diagonal of K. The importance of such scaling and why it helps has been explained for an 
idealized one-dimensional problem in [24]. We further reduce the condition number of the system matrix 
and the number of iterations for convergence by applying an incomplete Cholesky preconditioner with zero 
fill-in to the explicitly rescaled system, 

k w LL^. 

The finite element analysis in topology optimization requires the solution of a sequence of (usually) 
symmetric linear systems. In each optimization step, the algorithm updates the element densities, and after 
the first few optimization steps the changes in the design variables tend to be small from one optimization 
step to the next. Hence, the optimization leads to small changes from one linear system to the next, and the 
search space generated for one linear system provides important information for subsequent linear systems. 
First, the solution of one system can be used as an initial guess for the next system, reducing the initial 
residual. Second, an approximate invariant subspace derived from the Krylov space generated for one linear 
system can be used for subsequent linear systems, improving the convergence rate of the iterative solver. 
This is the basic idea of Krylov subspace recycling; however, other subspaces may also be used for 'recycling' 
[13) . Since the linear systems discussed in this paper are symmetric, we use the Recycling MINRES algorithm 
(RMINRES) for Krylov subspace recycling [53]. Unfortunately, solving a sequence of problems on meshes 
that change periodically makes recycling more difficult. Although it is not hard to map relatively smooth 
eigenvectors from a mesh to an updated mesh, the combination of mesh adaptation and preconditioning 
seems to give accuracy problems. Recycling is still effective for AMR; however, it is not nearly as beneficial 
as on a static mesh, and its improvement is work in progress. 

6 Implementation Issues 

For the implementation of adaptive mesh refinement, we use the libMesh library 11: developed at the 
University of Texas at Austin and the Technische Universitat Hamburg-Harburg. The libMesh library 
provides a C-I--I- framework for numerical simulations of partial differential equations on serial and parallel 
platforms. It supports one-dimensional, two-dimensional, and three-dimensional finite element and finite 
volume simulations on adaptive meshes. The libMesh software uses PETSc [2] [1] for the solution of linear 
systems on both serial and parallel platforms. However, we use our own custom linear solvers with Krylov 
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subspace recycling and preconditioners as detailed in Section [S] and references [HI [T3] ■ For compatibility 
with the libMesh package we have implemented the RMINRES method in the PETSc framework. For the 
incomplete Cholesky preconditioner we used routines provided by PETSc. 

We have developed two-dimensional and three-dimensional topology optimization algorithms on top of 
libMesh. Currently, we use element-based design variables, the SIMP method for material interpolation [S], 
the OC method for optimization [51 13], and Sigmund's filter technique [TBlfTOll?!] with some modifications for 
dealing with adaptive refinement. Following Stainko [22], we make a small modification in Sigmund's filter 
for a nonuniform mesh. Sigmund's filter takes a distance and density weighted average of the sensitivities of 
all elements in a certain radius as 



dc 
dpe 



1 



Pe l^d tide ^ Opd 



(3) 



where dc/ dpe is the sensitivity of the compliance with respect to the density of element e, and Hde is a 
distance weight defined as 

Hde = niax{r„,i„ — dist(c?, e), 0}. (4) 

The parameter r„i„ is a given radius for the filter (for the work reported here we use r„i„ = Tamr), and 
dist(d, e) is the distance between the centers of elements d and e. For a nonuniform mesh, we take the 
variation of element size into account by using the element volume to redefine the weight in the filter [22] as 



dc 

dpe 



PeJ2d tideVd 



Edc 
PdHdeVd-^ . 



(5) 



The filter radius r„,i„ is often a length scale independent of the mesh representation. Notice that the filter 
will be effectively deactivated if its size is smaller than that of the smallest element, i.e., no element has any 
neighbors within distance r„i„. 

Because of the hierarchical data structure of libMesh, we must allow level-one mesh incompatibility. 
However, we avoid level-two and higher mesh incompatibility. For example, if the configuration in Figure[3l[b) 
would result from mesh refinement, we refine the gray elements as well. If the configuration, in particular 
the gray elements, in Figure [3l[b) would result from mesh derefinement, we avoid the derefinement. In this 
way, we limit mesh incompatibility to level-one mesh incompatibility. As indicated by the circled nodes in 
Figure ini^a), level-one mesh incompatibility results in hanging nodes. The libMesh package handles those 
hanging nodes by using the projection method to enforce constraints in the stiffness matrix. We divide the 
degrees of freedom (DOFs) into two groups. Group one consists of all the unconstrained DOFs, and group 
two consists of the constrained DOFs on the hanging nodes. The constrained DOFs can be computed by 
linear interpolation from unconstrained DOFs. If we define vector u on the unconstrained DOFs, then 



u = 





I 







{pn) 


P 








(6) 



is the mapping of u to all the DOFs, where P is the interpolation matrix. We compute u by solving the 
projected system 

r r rtT 1 r T r, 1 r T rtT 1 

(7) 
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P 



P^ 




Since libMesh does not drop the constrained DOFs in the linear system, the projected system in ([7]) is 
singular when there is any hanging node. Krylov subspace methods can handle such singularities as long as 
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Figure 3: Mesh incompatibility with examples of quad, triangle and hex elements: (a) level-one mesh 
incompatibility marked by circled nodes; (b) level-two mesh incompatibility marked by circled nodes. We 
allow level-one mesh incompatibility (see Section [5]), but we avoid level- two and higher incompatibility by 
refining the gray coarse elements and by not derefining their children elements if these gray elements result 
from a potential derefinement. 



the right hand side is consistent, but these singularities may cause problems for preconditioners. To avoid 
the singularities, we set the diagonal entries in the matrix corresponding to the constrained DOFs to 1 and 
solve 
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In the end, we recover the constrained DOFs using the interpolation matrix: 



(8) 
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7 Results and Discussion 

We solve three problems to demonstrate the improvements of our new AMR scheme and verify that the 
computed designs using AMR meshes are equivalent to designs on uniform fine meshes. 

For the first (2D) test problem, we compute the optimal design on a uniform fine mesh and on an 
adaptively refined mesh with both our AMR scheme and an approach following references (HI [22]. The 
highest level of refinement in the AMR meshes has the same element size as that in the uniform mesh. 
The results show that our scheme computes a solution equivalent to the optimal design on the uniform fine 
mesh (within a small tolerance), while the alternative AMR approach from [H [12] does not. Moreover, the 
experiments elucidate how this suboptimal design arises from the strategy to only refine the results from a 
fixed coarse mesh. 

For the second (3D) test problem, we compare the optimal design using an adaptive mesh and our AMR 
strategy with the optimal design on a uniform fine mesh for a relatively simple cantilever beam problem. 
Again the maximum refinement level for the AMR mesh has the same element size as the uniform, fine mesh. 



10 



We also use this test to show that our AMR strategy leads to faster convergence in both the linear iterations 
(finite element solution) and the nonlinear iterations (design optimization) as well as to a significant reduction 
in runtime (about a factor of three). 

In the third test problem (also 3D), we compare the optimal design using an adaptive mesh and our 
AMR strategy with the optimal design on a uniform, fine mesh for a more complicated design problem. For 
all three test problems our AMR strategy leads to essentially the same design as is obtained on a uniform 
mesh, but at significantly reduced cost. 

To evaluate the relative difference between two designs, we need a quantitative measure. We define the 
the relative difference between two designs as 

We take p(l) to indicate the design on the uniform fine mesh, and p(2) to indicate the design on the AMR 
mesh. This difference can be computed in a number of ways. To simplify comparison, we take the uniform 
mesh to be the AMR mesh with maximum refinement at every point. So, we refine the AMR mesh to the 
same fine mesh level at every point (without changing the design), and then evaluate the 1-norm of the 
difference between the designs. 



Test 1: 2D cantilever beam 

We compute the optimal design for the 2D beam problem shown in Figure [U^a). We first compute the design 
on a uniform, fine mesh. Figure lUJb) shows an intermediate result, and Figure lljc) the converged design. 
Note how the truss member at the lower-right corner has risen up noticeably from the intermediate result 
to the final design. An effective AMR procedure must be able to capture such an evolution in the design. 

Next, we solve the same problem following the strategy mentioned in references [5J (35] • We start with 
a relatively coarse mesh (64 x 32), and obtain the converged solution to the topology optimization problem 
shown in Figure [Sfa). Then, we refine the mesh according to this coarse level result and we solve the 
optimization problem on this locally refined mesh until convergence, obtaining the solution shown in Figure 
[5ljb). Next, we refine the mesh and solve again. Finally, we obtain the result on the finest mesh shown in 
Figure [5]Jc). The truss member at the lower-right corner has remained roughly at its original position on the 
coarsest mesh in spite of the high resolution of the design, causing the resulting design to differ significantly 
from the optimal design obtained on the uniform, fine mesh, even though the smallest element size in the 
meshes in Figure [5jc) is the same as the element size for the uniform mesh shown in Figure S) This difference 
in material distribution is caused by the fine mesh discretization being confined to the region identified by 
the coarse mesh design. The mesh adaptation strategy only allows the fine mesh computation to refine the 
coarse mesh design. It does not allow the fine mesh computation to alter the design substantially, even if 
more accurate fine mesh computations indicate a better design. An additional problem is that the initial 
mesh needs to be relatively fine, such as the one in Figure [SJa), because a coarser initial mesh would lead 
to very poor solutions as the filter would be inactive at that mesh resolution. In this case, mesh adaptation 
with refinement only leaves fine elements that could be derefined for efficiency in the void regions. 

Now, we solve the same problem, starting with the coarse mesh of Figure [H^a) , but following our AMR 
strategy (Section |4]). We allow multiple mesh adaptations on any level and we maintain a layer of fine 
elements on the void side of the solid/void interface. This leads to the results shown in Figure [H with 
two intermediate results shown in Figures [HJa) and (b), and the final converged result in Figure |6l[c). Note 
how the truss member at the lower-right corner moves up as the optimization progresses, just as for the 
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Figure 4: Topology optimization on a 256 x 128 uniform mesh: (a) problem configuration (the volume 
constraint Vq is 50% of the domain volume); (b) an intermediate design; (c) final converged design. 
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(c) 

Figure 5: Topology optimization on an adaptive mesh with only refinement on each level: (a) converged 
result on the coarsest mesh with 2048 elements; (b) converged result on the intermediate mesh with 5675 
elements; (c) converged result on the final mesh with 20216 elements. Note the undesirable position of the 
truss member near the lower right corner, which remains nearly invariant during the evolution of the design. 
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(c) 

Figure 6: Topology optimization on an adaptive mesh with continual dynamic refinement and derefinement 
on each level: (a)-(b) intermediate designs; (c) final converged design on a nonuniform mesh with 25229 
elements, whose finest resolution is the same as the resolution of the uniform mesh in Figure |31 Notice that 
the truss member of the lower-right corner moves up as the AMR procedure progresses. 
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Table 1: Mesh adaptation scheme followmg Costa and Alves [8]. 



opt. step 


^max 


#eleni 


# unknowns 


relative Li difference 


80 


1 


2048 


4290 


21.16% 


125 


2 


5675 


11948 


19.53% 


200 


3 


20216 


41824 


19.42% 


300 


3 


20216 


41824 


19.60% 


400 


3 


20216 


41824 


19.66% 


500 


3 


20216 


41824 


19.66% 


600 


3 


20216 


41824 


19.65% 


700 


3 


20216 


41824 


19.67% 


800 


3 


20216 


41824 


19.66% 


900 


3 


20216 


41824 


19.65% 


1000 


3 


20216 


41824 


19.67% 



Table 2: Dynamic AMR scheme 



opt. step 


f 

^max 


#elem 


# unknowns 


relative Li difference 


27 


1 


2048 


4290 


21.13% 


37 


2 


6113 


12786 


19.80% 


100 


3 


24707 


50560 


17.77% 


200 


3 


25040 


51242 


13.00% 


300 


3 


25184 


51550 


5.00% 


367 


3 


25229 


51654 


0.17% 
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Figure 7: 3D cantilever beam example with domain scale 2:1:1. 

evolution of intermediate designs on the uniform mesh (reference). The figures also demonstrate how the 
mesh changes smoothly with the changes in material distribution. The smallest element size in the AMR 
meshes in Figure [6] is the same as the element size for the uniform mesh shown in Figure |4] (reference) . 
Compared with the final solution shown in Figure \E\ the solution obtained with our AMR strategy is closer 
to the solution obtained on the uniform mesh (reference). Indeed, based on the metric in (|10p . the relative 
difference between the designs in Figure [SKc) and Figure Sljc) is 19.6%, while the relative difference between 
the designs in Figure and Figure |4Kc) is only 0.168%. Furthermore, derefinement results in coarser 
elements in the void regions, cf. Figures [Sfc) and [SKc). These results are summarized in Tables [T] and 
where £niax refers to the highest refinement level present in the mesh. 

Test 2: 3D cantilever beam 

We compute the optimal design for the three-dimensional cantilever beam, shown in Figure[7l with a volume 
constraint of 25%. Exploiting symmetry, we discretize only a quarter of the domain. We solve this problem 
on a (fixed) uniform mesh with 128 x 32 x 32 B8 elements and also following our AMR strategy. The initial 
mesh for the AMR-based design has 64 x 16 x 16 B8 elements. The final results are shown in Figure |S] 
with Figure [HKa) displaying the solution on a uniform, fine mesh, and Figure [5{8) displaying the AMR 
solution; note the large blocks in parts of the void region in Figure |Hl[b) . The relative difference between 
these two designs is only 0.0909% (Eq. ([TU|). We use the preconditioned, recycling minimum residual solver 
(RMINRES) proposed in [24j to solve the linear systems arising from the finite element discretization for a 
given material distribution. The dimensions of the linear systems of equations for the adaptive mesh are less 
than half of those for the uniform, fine mesh. The difference is even larger early in the optimization iteration. 
Moreover, the number of RMINRES iterations for the linear systems derived from the adaptive mesh are 
slightly smaller than those for the uniform, fine mesh (much smaller early in the optimization), because 
the adaptive meshes tend to lead to better conditioned linear systems. Therefore, using AMR reduces the 
solution time roughly by a factor of three; see the statistics in Figure [9] 

Test 3: Cross-shaped domain 

We compute the optimal design for the more complex three-dimensional test problem shown in Figure 1101 
For the cross-shaped domain, we compute the optimal design subject to the fixed boundary on the bottom 
front and back ends, and two loads on the left and right sides. The maximum volume allowed is 20% of 
the domain volume. We solve this problem both on a uniform mesh and on an adaptive mesh following our 
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(a) 




(b) 

Figure 8: Final solutions of the 3D cantilever beam problem (Figure [T]) obtained using symmetry on a 
quarter of the domain as indicated by the mesh: (a) final solution on a fixed uniform mesh with 128 x 32 x 32 
elements; (b) AMR solution on a mesh with 57173 elements; the finest local resolution is the same as that 
of the uniform, fine mesh. 
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Figure 9: Comparison of linear solver statistics for the cantilever beam design problem on a uniform, fine 
mesh and on an adaptive mesh: (a) the number of unknowns in the linear systems arising from the finite 
element discretization; (b) the number of preconditioned MINRES and RMINRES(200,10) iterations (see 
below) for each optimization step; (c) solution times with MINRES and RMINRES(200,10) for the linear 
systems arising from finite element discretization at each optimization step. The parameters m and k 
in RMINRES(m,k) have the following meaning. The method recycles an approximate invariant subspace 
associated with the smallest k eigenvalues from one linear system to the next. In the solution of single linear 
system, the approximate invariant subspace is updated every m iterations. 
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Figure 10: A 3D compliance minimization problem in a cross-shaped domain with the bottom (shaded) front 
and back ends fixed and the bottom left and right ends pulled down. The volume constraint Vo is 20% of 
the domain volume. 

AMR strategy. The results are shown in Figures [Til and [T2| respectively. The uniform mesh consists of 40960 
B8 elements, while the final adaptive mesh consists of only 19736 B8 elements. Moreover, the optimization 
converges in over 200 steps on the uniform mesh, but in only 106 optimization steps on the adaptive mesh. 
The adaptive mesh refinement reduces the total solution time by more than a factor of three to about 30% 
of the solution time for the uniform mesh. Nonetheless, the relative difference between these two designs is 
only 2.58% (Eq. 1^). 

8 Conclusions 

In order to reduce the high computational cost of accurate three-dimensional designs by topology optimiza- 
tion we use adaptive mesh refinement. We propose several critical improvements to the approaches proposed 
by Costa and Alves jB] and Stainko [22j in order to attain better designs. In particular, we want to obtain 
the same optimal designs that would be obtained on a uniform, fine mesh with AMR discretization hav- 
ing significantly fewer elements but the same fine mesh resolution. The purpose of AMR is to reduce the 
cost for the (same) optimal design; we do not want to reduce the quality of designs. For large, complex, 
three-dimensional design problems we could not possibly use a uniform fine mesh at the desired resolution. 
Our approach requires a dynamic meshing strategy that involves continual refinements and derefinements 
following the strategy laid out in Section |31 Derefinements should also lead to further efficiency improve- 
ment by reducing the number of elements in void regions, especially for three-dimensional problems. Using 
three test problems, we demonstrate that our AMR algorithm achieves the desired designs that are within 
a small tolerance of those obtained on a uniform, fine mesh with the same finest level. Our AMR strategy 
significantly reduces the total runtime, nonlinear, and linear iterations with respect to using uniform meshes. 

Important future work includes error estimation in the finite element analysis and mesh refinement and 
derefinement governed by both considerations of accurate design and error estimation. In addition, we plan 
to work on preconditioners that can be adapted with the mesh (rather than recomputed) and to improve the 
convergence rate of Krylov methods with subspace recycling [211 [13] • We also intend to extend the present 
AMR technique to multiphysics problems ^Tj. 
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Figure 11: The optimal solution to the design problem shown in Figure [10] on a uniform finite element mesh 
with 40960 B8 elements. The quarter-mesh discretization is shown on the bottom figure. 
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Figure 12: The optimal solution to the design problem shown in Figure [TOl on an adaptively refined mesh. 
The final mesh consists of 19736 B8 elements. The quarter-mesh discretization is shown on the bottom 
figure. 
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